clc; clear;

%% 参数设置
alpha_list = linspace(2, 3, 100);    % alpha参数范围
q10_list = linspace(0, 1, 100);      % q10参数范围

A_template = [-1, 0, 0;
               0, 0, 0;
               0, 0, 0];
B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];

dt = 0.05;            % RK4步长
T_total = 600;        % 总积分时间
T_trans = 550;        % 暂态时间
N_steps = round(T_total/dt);
N_trans = round(T_trans/dt);

delta0 = 1e-8;        % 微扰大小

LE_vec = zeros(length(q10_list)*length(alpha_list),1);

fprintf('开始计算最大李雅普诺夫指数 (雅可比扰动轨道法)...\n');
parfor idx = 1:length(alpha_list)*length(q10_list)
    [i_q10, i_alpha] = ind2sub([length(q10_list), length(alpha_list)], idx);
    alpha = alpha_list(i_alpha);
    q10 = q10_list(i_q10);

    A = A_template;
    A(1,1) = -alpha;

    x0 = [1e-6; 0; 0; q10];

    LE = estimate_max_LE_Jacobian(@(x) mCNN_single(0,x,A,B), x0, dt, T_total, T_trans, delta0, N_steps, N_trans);

    LE_vec(idx) = LE;
end
fprintf('计算完成。\n');

LE_map = reshape(LE_vec, [length(q10_list), length(alpha_list)]);

%% 绘制热力图
figure;
h = imagesc(alpha_list, q10_list, LE_map);
set(gca,'YDir','normal');
colorbar;

% 设置NaN透明
alpha_data = ~isnan(LE_map);
set(h, 'AlphaData', alpha_data);

xlabel('\alpha');
ylabel('q_{10}');
title('最大李雅普诺夫指数 (雅可比扰动轨道法)');
grid on;

%% --- 基于雅可比扰动轨道法计算最大LE ---
function LE_max = estimate_max_LE_Jacobian(f, x0, dt, t_total, t_trans, delta0, N_steps, N_trans)
    h = dt;
    steps = N_steps;
    trans_steps = N_trans;
    
    x = x0;          % 主轨迹
    dim = length(x0);
    v = delta0 * [1; zeros(dim-1,1)];  % 固定扰动方向
    
    sum_ln = 0;
    count = 0;
    
    for i = 1:steps
        % 主轨迹RK4积分
        x = RK4(f, x, h);
        
        % 计算雅可比矩阵
        J = Jacobian(f, x, 1e-8);
        
        % 扰动向量线性演化
        v = v + h * (J * v);
        
        % 每10步归一化并累积log增长率
        if mod(i, 10) == 0
            dist = norm(v);
            if i > trans_steps
                sum_ln = sum_ln + log(dist/delta0);
                count = count + 1;
            end
            v = (v / dist) * delta0;
        end
    end
    
    LE_max = sum_ln / (count * h * 10);
    
    % 绝对值大于10设为NaN
    if abs(LE_max) > 10
        LE_max = NaN;
    end
end